function results = runEGARCHmodel(y)
cmaesOn     = 1;
params      = zeros(4,1);
params(1,1) = 0.01;    %intercept
params(2,1) = 0.80;     %persistence in log(sigma)
params(3,1) = 0.01;    %ARCH coefficient
params(4,1) = 0.01;    %ARCH coefficient
boundRoot   = 0.9999;

options     = optimset('Display','off','MaxIter',1e5,'MaxFunEvals',1e5,'TolFun',1e-12,'TolX',1e-12);
[paramsNM,fvalNM] = fminsearch(@likehood,params,options,y,boundRoot);

if cmaesOn == 1
    sigma              = 2;              %The step size
    insigma            = [0.001; 0.3; 0.3;0.2];
	%insigma            = [0.001; 0.3; 0.3;0.2;0.2];
    opts.SigmaMax      = 5;              %The maximal value for sigma
    opts.MaxIter       = 1e6;            %The maximum number of iterations
    opts.PopSize       = 50;             %The population size
    opts.VerboseModulo = 1D6;             %Display results after every 10'th iteration
    opts.TolFun        = 1e-8;           %Function tolerance
    opts.TolX          = 1e-8;           %Tolerance in the parameters
    opts.Plotting      = 'off';          %Dislpay plotting or not
    opts.Saving        = 'off';           %Saving results
    opts.Resume        = 'off';
    opts.Display       = 'off';    
    warning('off')
    [paramsCMAES,fvalCMAES] = cmaes_dsgeDisplay(@likehood,paramsNM,sigma,insigma,opts,y,boundRoot);
end
if fvalCMAES < fvalNM
    paramsOpt = paramsCMAES;
else
    paramsOpt = paramsNM;
end
    
% Computing standard errors
numParams = length(paramsOpt);
score     = zeros(length(y),numParams);
epsValue  = 1D-7;
for i=1:numParams
    params_peps = paramsOpt;
    params_peps(i,1) = params_peps(i,1) + epsValue;
    [~,~,logLike_peps]  = likehood(params_peps,y,1);
    
    params_meps = paramsOpt;
    params_meps(i,1) = params_meps(i,1) - epsValue;
    [~,~,logLike_meps]  = likehood(params_meps,y,1);

    score(:,i) = (logLike_peps-logLike_meps)/(2*epsValue);
end
A = zeros(numParams,numParams);
for k=1:length(y)
   A = A + score(k,:)'*score(k,:);
end
AVar = inv(A);

% Computing conditional vol. 
[sumLogL,logSigmaHat]  = likehood(paramsOpt,y,boundRoot);

% Computing an estimate of the variability in conditional vol, provided the
% sample is not too long (in which case it is simulated data)
rng(1,'twister')
numSim = 5000;
biasConVol2 = zeros(10,1);
if sum(~isnan(AVar(:))) && length(y)<500
    % Empirical distribution for the unconditional variability
    sig2Un      = mean(y.^2);
    AVar_sig2Un = (mean(y.^4)-sig2Un^2)/length(y);
    CI_sig2Un = [sig2Un-5*sqrt(AVar_sig2Un),sig2Un+5*sqrt(AVar_sig2Un)];
    
    AVarChol  = chol(AVar,'lower');
    sigmaHatSim = zeros(length(y),numSim);
    sumConVol2  = zeros(10,numSim);
    s           = 1;
    while s <=numSim
        paramsMC = paramsOpt + AVarChol*randn(numParams,1);
        [sumLogLMC,logSigmaHatMC]  = likehood(paramsMC,y,1);
        % checking for extreme outliers inconsistent with unconditional distribution for y.^2
        sig2UnMC = mean(exp(2*logSigmaHatMC));
        if isnan(sumLogLMC) || sig2UnMC>CI_sig2Un(1,2) || sig2UnMC<CI_sig2Un(1,1)
            
        else
            sigmaHatSim(:,s) = exp(logSigmaHatMC);
            sumConVol2(1,s) = sigmaHatSim(:,s)'*sigmaHatSim(:,s);
            for j=1:9
                sumConVol2(1+j,s) = sigmaHatSim(1:end-j,s)'*sigmaHatSim(1+j:end,s);
            end
            s = s+1;
        end
    end
    biasConVol2(1,1) = mean(sumConVol2(1,:),2) - exp(logSigmaHat)'*exp(logSigmaHat);
    for j=1:9
        biasConVol2(1+j,1) = mean(sumConVol2(1+j,:),2) ...
            - exp(logSigmaHat(1:end-j,1))'*exp(logSigmaHat(1+j:end,1));
    end
else
    sigmaHatSim = NaN(length(y),numSim);    
end
% Saving output
results.params       = paramsOpt;
results.paramsSE     = sqrt(diag(AVar));
results.params_tstat = paramsOpt./sqrt(diag(AVar));
results.conVolSim    = sigmaHatSim;
results.conVol       = exp(logSigmaHat);
results.biasConVol2  = biasConVol2;
results.sumLogL      = -sumLogL;
results.AVar         = AVar;
warning('on')


    function [sumLogL,logSigma,logLike] = likehood(params,y,boundRoot)
        % test for stationarity
        if abs(params(2,1))>=boundRoot
            sumLogL = NaN;
            logSigma = NaN;
            logLike = NaN;
            return;
        end
        T            = size(y,1);
        logSigma     = zeros(T,1);
        residuals    = y(1:T,1);
        logSigma(1,1)= params(1,1)/(1-params(2,1));
        logLike      = zeros(T,1);
        for t=1:T
            logLike(t,1)  =  - logSigma(t,1)-0.5*(residuals(t,1)/exp(logSigma(t,1)))^2;
            if t < T
                logSigma(t+1,1) = params(1,1) + params(2,1)*logSigma(t,1)...
                    +params(3,1)*residuals(t,1)/min(exp(logSigma(t,1)),exp(-8))...
                    +params(4,1)*(abs(residuals(t,1)/min(exp(logSigma(t,1)),exp(-8))) -sqrt(2/pi));
                % Note E[abs(z)] = sqrt(2/pi)
            end
        end
        % We use a minimization routine
        sumLogL = -sum(logLike);
    end

end

